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Abstract. - We address the problem of understanding the variable abundance of 3- node and 
4-node subgraphs (motifs) in complex networks from a dynamical point of view. As a criterion 
in the determination of the functional significance of a n-node subgraph, we propose an analytic 
method to measure the stability of the synchronous state (SSS) the subgraph displays. We show 
that, for undirected graphs, the SSS is correlated with the relative abundance, while in directed 
graphs the correlation exists only for some specific motifs. 



Recent empirical evidences indicate that complex networks, among other common prop- 
erties, are characterized by the presence of various length cycles and specific motifs [1-3]. 
A motif M is a pattern of interconnections occurring either in a undirected or in a directed 
graph G at a number significantly higher than in randomized versions of the graph, i.e. in 
graphs with the same number of nodes and links (and eventually degree distribution) as the 
original one, but where the links are distributed at random. As a pattern of interconnections, 
is usually meant a small connected (undirected or directed) n-node graph M which is a sub- 
graph of G. The concept of motifs was originally introduced by U. Alon and coworkers, who 
studied motifs in biological and non-biological networks [4-8] . The research of the significant 
motifs in a graph G is based on matching algorithms counting the total number of occurrences 
of each n-node subgraph M in the original graph and in the randomized ones. The statistical 

significance of M is then described by the Z -score, defined as: Zm = M$ — ) where #m 

#M 

is the number of times the subgraph M appears in G, and ) and o#^J are, respectively, 

the mean and standard deviation of the number of appearances in the randomized network en- 
semble [5,7]. The reasons of the variable frequency of different n-node subgraphs in a specific 
network are still poorly understood. There are at least two possible explanations. On the one 
hand, it is possile that certain constraints on the growth mechanism of a network as a whole 
determine which motifs become abundant [9,10]. On the other hand, it is well known that the 
structure has important consequences on the network dynamics and functional robustness. So 

© EDP Sciences 



2 



EUROPHYSICS LETTERS 



that a particular n-node graph can become overrepresented because, due to its structure, it 
possesses some relevant functional properties [5] . 

In this letter, we address the question of network motifs in biological networks from a 
dynamic systems point of view. Naturally, a comprehensive analysis of the dynamics of 
networks is considerably more complicated than the corresponding analysis of their structure. 
This is due to the potentially complex functional dependencies between nodes, and to lack 
of knowledge of the specific interaction parameters. For such a reason, instead of modeling 
in details one particular biological network, we analyze the generic dynamic properties that 
arise from the topology of a n-node graph. In particular, we focus on the emergence of 
collective dynamic behaviors, such as synchronization, that is relevant in many biological 
systems, and we propose an analytic method to estimate the stability of the synchronous 
state (SSS) displayed by a n-node graph. We finally show that the SSS, potentially, can help 
explaining why certain network motifs are overrepresented in some real biological networks, 
while others are not. 

We assume that the dynamics of a n-node motif M can be represented as a system of n 
ODE's: 

Xi = fj(xi,X2,...,x n ) i = l,...,n (1) 

where xj G R m is the ?n-dimensional vector describing the state of node i (for instance the 
concentration of molecule i in a metabolic reaction, or the polarization state of neuron » in a 
neural network), and : R mxn — > R m is the function representing the effects on Xj of all the 
nodes connected to i. In particular, we are neglecting the influence of the other nodes of the 
graph G on the n-node motif M, and we are assuming that the fj's do not contain an explicit 
dependence on time. The issue of the stability of the steady states of Eqs. ((TJ), i.e. the sets of 
values (x*,X2, . . . ,x*) such that xj = x^ = . . . = x* = 0, has been investigated in Ref. [11]. 
Here we focus on the stability of the synchronized dynamics of Eqs. Jl}, which can be treated 
analytically within the context of the so-called Master Stability Function approach [3, 12-14]. 
In particular, we restrict to the case in which the equations of motion can be written as: 

x, = Fj(xi) + crJ2™ =1 Oii[Hjj(xj) - Hii(xi)] i = l,...,n. (2) 

where F, (x) : R m — > R m is the function governing the local dynamics of node i, Hij(xj) : 
R m — > R m describes the influence of node j on node i, a > is the coupling strength, and 
dij are the elements of the n x n adjacency matrix of graph M. In the case of a undirected 
n-node motif M, a^j = aji = 1 iff there is an edge joining node i and node j, and a^j = aji = 
otherwise. In the case of a directed motif, we assume o,j = 1 iff there is a directed edge from 
node j to node i, while = otherwise. Equations ^ can be rewritten as: 

ii = F^Xj) - a YTj=i ! >j H y( x j) i = l,...,n. (3) 

where £y = 5»j (XZz a n) ~ a ij are ^ ne elements of a zero row-sum (y]j hj — Vi) n x n matrix 
L with strictly positive diagonal terms (la > Vi). In the case of a undirected motif M, L is 
symmetric and coincides with the standard Laplacian matrix of the graph M [3] . In the case 
of a directed graph, the off-diagonal elements hj of L are respectively equal to — aij, while 
the i-th diagonal entry is equal to the in-degree of node i, fc| n = J2i a u- ^ n or der to proceed 
with the analytic treatment, we make the explicit assumption that the network is made of n 
identical and identically coupled dynamical systems. This corresponds to take in Eqs. ([2]) and 
Eqs. (J2J) Fj(xj) = F(x) Vi, and Hy (xj) = H(x) Vi,j. This assumption and the fact that L 
is zero-row sum, ensure the existence of an invariant set X! (t) = X2 (t) = • • • = x„ (t) = x s (t) , 
representing the complete synchronization manifold S. The main idea, first proposed by 
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Pecora and Carrol [12], is that the linear stability analysis of the synchronized state of Eqs. © 
can be divided into a topological and a dynamical part [12, 14]. Since the coupling term of 
Eqs. ([3]) vanishes exactly on S, a necessary condition for the stability of the synchronous state 
is that the set of [n — 1) * m Lyapunov exponents corresponding to phase space directions 
transverse to the synchronization manifold are entirely made of negative values. Considering, 
then, the mxn column vectors X = (xi, X2, . . . , x„) T and <5X = (<5xi, . . . , Sx n ) T , where 
Sxi (t) = Xi (t) — x s (t) is the deviation of the i^ h vector state from the synchronization manifold, 
one gets the variational equation: 

SX = [I n ® JF(x s ) -crL® JH(x s )]5X, (4) 

where I n is the nx n identity matrix, ® stands for the direct product between matrices, and J 
denotes the Jacobian operator. The first term in Eq. is block diagonal with mxm blocks, 
while the second term can be treated by diagonalizing L. 

We first concentrate on the case of undirected motifs, i.e. on symmetric and thus diag- 
onalizable laplacian L. Let Xi be the set of n real eigenvalues of L (Lvi — AjVj, i = 1, . . . , n), 
and Vj the associated orthonormal eigenvectors (vj • = 5ij). If L is symmetric, all its eigen- 
values are real, and they can be ordered by size as: = A x < A 2 < . . . < X n . The arbitrary 
state SX. can be written as 5X = J27=i v * ® C*(*)i where Q = (&,», ■ ■ ■ , Qm.i)- By substituting 
into Eq. (HJ, and using the condition that the eigenvectors are linearly independent, one is 
finally left with a block diagonalized variational equation, with each of the n blocks having 
the form of a variational equation for the coefficient (k(t)'- 

^-=K k ( k , k = l,...,n (5) 

where Kj, = [JF(x s ) — o-AfeJH(x s )] is the evolution kernel. Each equation in © corresponds 
to a set of to conditional Lyapunov exponents along the eigenmode corresponding to the 
specific eigenvalue Afc. For fc = 1, Ai = 0, and we have the variational equation for the 
synchronized manifold S. The m corresponding conditional Lyapunov exponents equal those 
of the single uncoupled system x = F(x), therefore no conditions on them will be imposed 
(in principle, the synchronized state itself can well have positive Lyapunov exponents and be 
chaotic). 

Notice that the Jacobian JF(x s ) and JH(x s ) are the same for each block k, since they are 
evaluated on the synchronized state. Consequently, the form of each of the blocks in Eqs. ([5|) 
is the same, with the only difference being in the multiplier A^. This leads one to replace aXk 
by v in Eq. ((5|), and to consider the generic m-dimensional variational equation: 

C = K,C = [JF(x s ) - !/JH(x s )] C, (6) 

from which one can extract the set of m conditional Lyapunov exponents as a function of the 
real parameter v > 0. The parametrical behavior of the largest of such exponents, A(^), is 
called Master Stability Function [12-14]. In fact, given a coupling strength er, one can locate 
the point aXk on the positive v axis, and the sign of A at that point will reveal the stability 
of that eigenmode. If A(crAfc) < Vfc = 2, ...,n, then the synchronous state is stable at the 
coupling strength a. 

In order to evaluate whether the stability of the synchronous state is favoured by the 
topology in a given 77-node graph more than in another, we adopt the following measures of 
stability. First, we assume that K(v = 0) > 0, meaning that the uncoupled systems x = F(x) 
support a chaotic dynamics. For v > 0, there are three possible behaviors of A(i/), defining 
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throe possible classes for the choice of the functions F(x) and H(x). Case I (II) corresponds 
to a monotonically increasing (decreasing) A(z/). Case III admits negative values of A(v) in 
the range v Cl < v < v C2 (see Fig 5.1 of Ref. [3]). For systems in class I, one can never stabilizes 
synchronization in any graph topology. In fact, for any a and any eigenvalues' distributions, 
the product a\k always leads to a positive maximum Lyapunov exponent, and therefore the 
synchronization manifold S is always transversally unstable. Class II systems always admits 
synchronization for a large enough a. In fact, given any eigenvalue distributions (any graph 
topology) it is sufficient to select a > v c /\2 (A2 7^ in a connected graph [15]) to warrant that 
all transverse directions to S have associated negative Lyapunov exponents. The synchronous 
state will be stable for smaller values of a in a graph with a larger A2, so that A2 can be 
used as a measure of the stability of the synchronous state (SSS). For systems in class III, 
the stability condition is satisfied when a > v Cl j\i and a < v C2 /Xn, indicating that the more 
packed the eigenvalues of L are, the higher is the chance of having all Lyapunov exponents 
into the stability range [14]. Consequentely, the ratio A2/A„ can be used as a measure of 
SSS. Classes II and III include a large number of functions F, describing several relevant 
dynamical systems, as the Lorenz and Rossler chaotic oscillators, and the Chua oscillator. It 
is important to notice that not only F, but also H has a role in determining to which class 
a specific dynamical system belongs to. As an example, a nearest neighbor diffusive coupling 
on the Rossler chaotic system yields a class II (class III) Master Stability Function, when the 
function H extracts the second (the first) component of the vector field [16]. In Fig. Q] (panel 
a and b) we report the two indices of SSS, namely A2 (class II) and A2/A4 (class III), for the 
six 4-node undirected motifs. We observe a general increase in the SSS's as the number of the 
edges in the motif increases. Such an increase in SSS is in agreement with the decrease of the 
synchronization threshold observed numerically in the Kuramoto model by Moreno et al. [17]. 
The two measures of SSS we propose are also in good agreement with the natural conservation 
ratio (NCR) for the same 4-node motifs in the the yeast protein interaction network reported 
in panel c). The NCR is a measure proposed in Ref. [18] to quantity the conservation of a 
given motif in the evolution across species, and is highly correlated to the motif Z-score. In 
panel d) and e) we show that SSS's and NCR are linearly correlated: we have obtained a 
correlation coefficient respectively equal to 0.94 and 0.93. This is an indication that motifs 
displaying an improved stability of cooperative activities (as synchronous states) are preserved 
across evolution with a higher probability. 

We now turn our attention to directed motifs. In a directed graph, the matrix L is asym- 
metric and in general not always diagonalizable. Nevertheless, L can be transformed into a 
Jordan canonical form, and it has been proven that the same condition valid for diagonal- 
izable networks (A(cAfc) < Vfe = 2, ...,n) also applies to non-diagonalizable networks [19]. 
In addition, the spectrum of L is either real or made of pairs of complex conjugates. Be- 
cause of the zero row-sum condition, L always admits Ai = 0, and the other eigenvalues 
Afe = A^ + iA^, k = 2, . . . ,n (having non negative real parts according to the Gerschgorin's 
circle theorem [20]) can be ordered by increasing real part (0 < X§ < ■ ■ ■ < A^). Conse- 
quently, the parametric equation (|6|) has to be studied for complex values of the parameter 
v = v + iv 1 . This yields a master stability function A(i/) as a surface over the complex 
plane i>, that generalizes the plots for the case v real. By calling 1Z the region in the com- 
plex plane where A(^) provides a negative Lyapunov exponent, the stability condition for the 
synchronous state is that the set {a\k,k = 2, . .. , n} be entirely contained in 1Z for a given 
a. This is best accomplished for connection topologies that make Af 1 as large as possible 

for class I systems, and for topologies that simultaneously make y|- as large as possible and 
max{| \\, |} as small as possible, for class II systems. 
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Fig. 1 - The value of SSS for each of the six 4-node undirected motifs are reported in panel a) for class 
II systems and in panel b) for class III, and compared with the natural conservation rates (NCR) in 
the yeast protein interaction network [18]. The motif identification number is the same as in Ref. [7]. 
In panel d) and e) we plot the values of SSS as function of the NCR (symbols) , and the linear fittings 
obtained (dashed lines). 

In Fig. [21 we consider the thirteen 3-node directed motifs. Two of them, namely motifs 
#3 and motif #11 give rise to non-diagonalizable L. Motif #8 is the only case where the 
eigenvalues are not real. In the left (right) panels we report for class II systems (Af /A^ 
for class III systems). The SSS measures are compared with the Z-score profile obtained for 
five different real biological networks, and shown as hystograms in the figure. Both class I 
and class II systems exhibit an average increase of SSS as a function of the number of links in 
the motif. However, the overall agreement of the SSS and the Z-score profiles is not as good 
as in the case of undirected 4-motifs. Here, we have obtained rather small values (ranging 
from 0.1 to 0.3) of the correlation coefficient, with a better agreement found in the case of 
the STKE network (panels c), Drosophila (panels d) and C.elegans (panels e), rather than 
in the transcriptional regulatory networks (panels a and b). This might be due to the fact 
that synchronization processes are more important in neural systems than in other biological 
systems as transcriptional networks, especially the simplest ones (E. coli and S. cerevisiae). 
We have also reported in figure, as dashed lines, the measure of the stability of stationary 
states proposed by Prill et al. [11]. Such a measure seems to be better indicated for those 
systems where the stability of stationary states can be a more relevant dynamical quantity 
to investigate than the stability of synchronous states. Fig. [2] clearly indicates that in some 
motifs, SSS and Z-score are better correlated than in others. Hence, for each motif M, we 
have defined an overlap coefficient Om as: Om = SSSm x %m- The maximum possible value 
Om = 1 indicates a perfect correlation between SSS and Z-score. The overlap coefficients 
obtained for the five studied systems are reported in Fig. [3] For both class II and class III 
systems we have high values of the overlap for motifs: 1, 7, 10, 12. 

Finally, we have considered the 199 4-node directed motifs. Here we report the results 
for three of the most statistically relevant motifs found in biological networks: the bifan, the 
biparallel and the feedback loop (see Ref. [5]). Such three motifs correspond all to cases in 
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Fig. 2 - The SSS of each of the thirteen 3-node directed motifs is reported (continuous line) for class 
II (left panels), and class III systems (right panels). For class III system the SSS values has been 
normalized to the maximum value so to vary in the range [0,1]. The SSS values are compared with 
the Z-score (hystograms) and with a measure of the stability of stationary states (dashed line) from 
Ref. [11], in five different biological networks: the transcriptional regulatory networks of E. coli (panels 
a) and S. cerevisiae (panel b), the signal transduction knowledge environment (STKE) network (panel 
c), the developmental transcriptional network of Drosophila melanogaster (panel d), and the neural 
connection map of C. elegans (panel e). 



which L can be diagonalized. The biparallel graph, that is abundant in the C. elegans and in 
transcriptional networks, has real eigenvalues and a relatively high value of SSSs: A^ = 1 and 
Af/A^ = 0.5. The same is true for the 4-node feedback loop (also found abundant in electric 
circuits [5]), having \f = 1, Af/A^ = 0.5 and maxfc> 2 {| Xf. |} = 1. Conversely, the bifan is 
not compatible with synchronization for any choice of F(x) and H(x), and for any value of 
cr, since A2 = and we have assumed the case of networked chaotic systems (A(u = 0) > 0). 
In fact, A^ 7^ iff the graph embeds an oriented spanning tree, (i.e., there is a node from 
which all other nodes can be reached by following directed links) [19,21] and this condition, 
that generalizes the notion of connectedness for undirected graphs [15] to directed graphs, is 
not valid in the case of the bifan. 

We warmly thank R.J. Prill and A. Levchenko for having provided us with their results 
on the stability of stationary states, and G. Russo for useful comments. S.B. acknowledges 
the Yeshaya Horowitz Association through the Center for Complexity Science. 
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Fig. 3 - The overlap coefficients for each of the thirteen 3-node directed motifs, and the five biological 
networks considered, are reported for class II (panel a) and class III systems (panel b). 
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